Orthonormal Bernstein Galerkin technique for computations of higher order eigenvalue problems

The numerical approximation of eigenvalues of higher even order boundary value problems has sparked a lot of interest in recent years. However, it is always difficult to deal with higher-order BVPs because of the presence of boundary conditions. The objective of this work is to investigate a few higher order eigenvalue (Rayleigh numbers) problems utilizing the method of Galerkin weighted residual (MWR) and the effect of solution due to direct implementation of polynomial bases. The proposed method develops a precise matrix formulation for the eighth order eigenvalue and linear electro-hydrodynamic (EHD) stability problems.• The article explores the same for tenth and twelfth order eigenvalue problems.• This method involves computing numerical eigenvalues using Bernstein polynomials as the basis functions.• The novel weighted residual Galerkin technique's performance is numerically validated by comparing it to other numerical/analytical approaches in the literature.


Introduction
Higher even order BVPs occur in various physical and biological problems. Hydrodynamics, hydro magnetic stability and astrophysics are a few to name. Many researchers have investigated higher even order BVPs due to their scientific significance and applications in a variety of fields of applible branches of sciences. Based on the literature review, many researchers have attempted to address higher order boundary value problems. Nevertheless, only a few studies have investigated higher order differential eigenvalue problems that arise in hydrodynamic and hydro-magnetic stability theory. Eighth, tenth, and twelfth order eigenvalues have endured extensive numerical computation in [ 1 , 2 ]. Additionally, the authors [3] have explored an eigenvalue problem in EHD that indicates the presence of electric forces and is modeled by an 8th order differential operator. The explicit use of numerical methods in some of these scenarios may produce erroneous results because of the bifurcation issues with the stable solutions of well-known Navier-Stokes equations (or of certain existing models). It is to be noted that bifurcations occur as the eigenvalues depend on various physical parameters. Consequently, an analytical solutions as well as numerical investigation of hydrodynamic stability theory is greatly desired.
When the bottom layer of fluid is highly agitated underside under the influence of rotation, the numerical computation of eighth, tenth, and twelfth order [1][2][3] eigenvalues occur. The top layered fluid will be heavier than that of the fluid at the bottom and in this state, which causes the layers to be potentially unstable. Thus, viscosity on the part of the fluid has a prevent tendency to rearrange itself. In this work, we are influenced by any supplementary effects of spins, and their rotation will initiate new components into the resulting thermal instability. In fluid flow, instability becomes entrenched typically as over stability in the presence of rotation, but it continues as stationary convection due to the impact of a magnetic field [ 1 , 2 ]. When the instability of a magnetic field manifests as over stability, eighth order differential eigenvalue difficulties arise [1] and can be illustrated in the following form: The Rayleigh critical value, which is the smallest eigenvalue, changes as gravity changes [2] . In Straughan's monograph [4] , the hydrodynamic stability equations are described in further detail. The lowest possible value of and attained by solving Eqs. (1) and ( 1a ) for the matched value of , which may be complex, it is supposed that , is complex. However, for the Engineering significance of , compels it to be real and subsequently ( = ) is purely imaginary [1] .
When a homogeneous magnetic field over the fluid in the same direction as gravity induces ordinary convection and over stability, respectively, tenth and twelfth order differential eigenvalue issues [1,2] arise. The essential condition for stability is popularly known as Rayleigh criterion. Tenth order differential eigenvalue problem which occurs as instability is of the form [1] .
Twelfth order equation takes the form [1] ( The stability of fluid flow is determined by numerical value of the non-dimensional parameters, referred to as Taylor numbers which gives a measure of extent to which Rayleigh's criteria is violated.
This study utilizes a set of orthonormal Bernstein polynomials over a finite domain. It is tempting to implement due to the convergence of Bernstein polynomials for specific features, as well as non-negative optimal stable bases among other bases [5][6][7][8][9] . Considering functions of bounded variation, the authors of the study [5] estimated the convergence rate of Bernstein polynomials in connection with the arithmetic means through the series of overall variation. Several researchers have used Bernstein polynomials in recent years as a practical tool because of their optimal stable bases among nonnegative bases and their dependability for a variety of noteworthy features used in computer-aided geometric design [5] . The authors [6] employed operational matrices of Bernstein polynomials and its expansion to solve parabolic PDEs. The studies [ 7 , 8 ], provide numerous intuitive geometrical qualities, elegant algorithms, and great numerical stability under the veil of Bézier curve and surface representations [9] gives an overview of Bernstein polynomial applications and illustrates how the Galerkin method is used to solve high even-order differential equations.
In general, an even order Benard type convection problems (the 2 k -th order, k = 4, 5, 6) together with homogeneous BCs can be written as [9] ( 2 ) + If the stated boundary conditions below are met There are a few works of literature on the approximate eigenvalues of higher order Sturm-Liouville problems. The authors [1] focused on two numerical approaches of finite difference method, namely direct numerical technique, second and fourth-order finite difference techniques to compute eigenvalue problems illustrated therein. In [1] , they demonstrated that their direct numerical approach (without reducing the higher order derivatives) eventually brings an ill-conditioned system and undergoes word length difficulties as well. The authors updated the various existing software to compute the eigenvalues with good accuracy. For the specified values of the parameters ( , , ) , the computed each eigenvalue becomes complex. The authors varied a and to minimize the imaginary part. Besides, direct methods cause some difficulties due to software implementations. As a result, some infinite eigenvalues occur. However, their second approach avoids this difficulty by giving complex banded coefficient matrices, but much more expensive to implement. From these discussions, we aim to develop a simple and efficient technique which works well for these kinds of problems with less difficulties. The thermal instability of fluid flow arises in sixth order Benard convection problems and Rayleigh numbers have been studied thoroughly in [10][11][12] applying numerous spectral techniques.
Several numerical methods, namely local adaptive differential quadrature have been implemented in various studies. A numerical method namely Adomian decomposition methods [ 13 , 14 ], energy inequalities method [15] , Galerkin based Septic B-Splines technique [16] , splines and non-polynomial spline techniques [17][18][19] , Local adaptive differential quadrature method [20] , reproducing kernel Hilbert space method [21] to solve the boundary value problems of higher order with multi-boundary conditions. The authors of [3] used both analytical and numerical methods to investigate the linear Electro-hydrodynamic stability problem of an eighth order differential equation. The latter one is implemented through spectral Galerkin and collocation methods taking Chebyshev and Legendre polynomials [12] as basis functions. The rigorous study of eighth, tenth, and twelfth order linear and nonlinear BVPs using Galerkin MWR and polynomials as the basis function is found in [22][23][24][25] . The numerical approximations on eighth, tenth and twelfth order eigenvalue problem is infrequent. Therefore, we focused on an efficient numerical approach for solving equations analogous to (4). A general approach to solve an even-order Sturm-Liouville problem has recently been developed [26] using the Lie group method and the Magnus expansion (tested for up to eight). Thus to the best of our knowledge a general higher even order such problems has not been studied before using Orthonormal Bernstein Galerkin Technique. This fact motivates us to study further focusing Eq. (4) using Bernstein Galerkin Technique to test the efficiency of considering such polynomials.
In the proposed method, basis functions are satisfied by the Dirichlet-type boundary conditions. On the other hand, in the weak form of the Galerkin residual equation, all the essential types of boundary conditions are incorporated directly. The vanishing at the two boundaries of the interval gives greater flexibility to the Bernstein polynomials and is found to be more attractive for implementing the Galerkin method of weighted residual. Difference schemes are used in [27] to solve the one-dimensional wave equation with nonlocal boundary conditions. The circuit's time and frequency domain properties are examined [28] , the filter parameters are chosen analytically, and the outcomes are numerically tested.
The rest of the article is organized as follows: Several well-known properties of Bernstein polynomials shifted orthonormal Bernstein polynomials are listed in section 2. Section 3 illustrates the Bernstein Galerkin method to solve higher even order eigenvalue problems. Section 3 further explores error estimation and the convergence of the residual function. To confirm the reliability and proficiency of the presented method, several numerical examples of BVPs, available in the literature, are presented in Section 4. Section 5 of this study mentions its findings as a short conclusion.

Bernstein polynomials
Even though Bernstein polynomials have numerous advantageous characteristics, they lack the attribute of orthogonality. The orthogonality property is particularly helpful for many applications, including least squares approximation and finite element methods. As a result, Bernstein polynomials are frequently less practical to use in these methods than conventional orthogonal polynomials like Legendre, Chebyshev, or Jacobi polynomials.
The optimal stable, nonnegative Bernstein polynomial basis gained popularity among several researchers and useful properties have been applied to solve regular as well as singular second to higher order BVPs and eigenvalue problems [ 9-11 , 22 , 23 , 29-38 ]. As stated in [22] , the polynomials can be used to form a basis over [0, 1] which is complete. Following are the specifications for Bernstein polynomials over the finite range [0, 1] where n being positive integer and

A few important and useful properties of Bernstein polynomials
All Bernstein polynomials are zero at the end points in [0, 1], that is, The n -th derivative of Bernstein basis can be stated as (vi) The binomial expansion of the nth degree Bernstein polynomial can be simply expressed in terms of any higher power base

Shifted orthonormal Bernstein polynomials
Orthonormal Bernstein polynomials over [0, 1] can be explicitly represented by utilizing the Gram-Schmidt procedure on the set of Bernstein polynomials in Eq. (3) of different degree n is given as Utilizing Bernstein non-orthonormal basis, the above expression takes the form The orthonormal Bernstein polynomial , ( ) is the th eigenfunction of the singular Sturm Liouville problems In terms of the weight function, these polynomials have an orthogonality relationship i.e., ( ) = 1 .
The orthonormal Bernstein polynomials necessarily satisfy the following relationships over the interval [ 0 , 1 ]

Orthonormal Bernstein-Galerkin technique
Here we move on to the main goal and focus on our key issues in this section. To be specific we are now interested in finding the solution to Eq. (4) that also satisfies boundary condition 4(a) using Bernstein polynomials. The integrals to the higher order eigenvalue problem defined in (4) is considered here.
The normalized Bernstein Galerkin approximation [9] to the Eq. (4) is to obtain ∈ in a way such that, where the inner product ⟨ , ⟩ = ∫ ( ) ( ) ( ) is defined as in 2 ( ) and its norm It is of utmost importance to highlight that choosing a suitable basis for to achieve the Bernstein-Galerkin approximation to (10) that yields the simplest linear system is the key challenge in using the Galerkin Bernstein method of weighted residual. An appropriate normalized basis for Galerkin MWR can be written as where ( ) for all = , + 1 , ....., − . Because of the 2k boundary conditions, the first and last k expansion coefficients are both zero.

Convergence of Galerkin MWR
Stability and convergence measures [9] for the polynomial estimations in the Galerkin MWR have been studied by several authors [ 9 , 10 , 22-25 , 31-33 ]. Moreover, orthogonal polynomials basis can approximate general functions in Hilbert spaces. Here, we reveal a short appraisal on this topic.
Function approximation: Here we recall Bessel inequality ( Theorem: (Bessel inequality) Let ( ) be an orthonormal sequence in an inner product space i.e., in H , then for every ∈ H  9) . We can define as Theorem. [39] : Let be an arbitrary element in . Since is a finite dimensional complete and closed subspace, therefore, is a subset of . Therefore, has the best unique approximation out of , say Since, ∈ , there exists a unique coefficient ] .
Let the exact and the approximate solutions be ( ) and ̃ ( ) , respectively, provided the functions ( ) are linearly independent. For detailed demonstration on convergence, we refer our previous work [ 10 , 11 ] The sequence of estimated eigen-pairs converge to the precise integrals in (20) as the number of degrees of freedom grows indefinitely. This compels an uniform rate of convergence at every point in the domain. The eigen solution Eq. (19) must be used to optimize (minimize) the residual error Now we move on to demonstrate a theorem for the convergence of the suggested approach.

Numerical illustrations
To assess the competency of the suggested method, four linear eigenvalue problems with two classes of boundary conditions have been worked out. The eigenvalues computed by the current technique are compared with the various numerical and analytical techniques [ 1 , 2 , 14 , 26 ] available in the literature. The Galerkin method is tested with finite intervals from eight to twelfth orders along with regular endpoint boundary conditions. Example 1. We study the eighth order Sturm-Liouville problem [26] ( 1 ) = ′′ ( 1 ) = ( 4 ) ( 1 ) = ( 6 ) with exact eigenvalues 2 = ( ) 8 The relative errors of the first five eigenvalues using 15 Bernstein polynomials ( n = 15) converge and match with the Magnus method [26] in the interval shown in Table 1 . If the number of polynomials is increased even further, stable convergence is maintained. Our current method can accurately determine the first eigenvalues of the problem under investigation.

Example 2.
We investigate the linear stability of the steady state response in an electrohydrodynamic convection model in a layer situated toward the normal mode of perturbations [1] as defined by [3] and placed between the walls at = ±0 . 5 .
The boundary conditions for derivatives of even order are illustrated as The general solution can be found using the direct method by determining the determinant of the characteristic equation and being subjected to the multiplicity of the roots of the equation relevant to the problem [3] . Here R refers to the Rayleigh number, and satisfy the "ellipticity " condition given by 4 The authors [3] applied direct analytical/numerical schemes and thus completed the investigation of this problem. It is evident from their effort that the Rayleigh number occurs as the lowest one which is real and positive i.e., > 0 , and satisfies the eigenvalue problem (26a). If a, L , ≠ 0 , the discussions of a multiplicity of the eigenvalues become difficult. So, the use of the direct analytical method is quite incomprehensible, and an alternative numerical method is sought.
Following the above facts, the authors' in [3] applies analytical method and Chebychev and shifted Legendre polynomials spectral methods to formulate the Eq. (26a) explicitly. The authors' used 2 strategy worked out in their previous article [12] .
The numerical evaluations of the first nine critical Rayleigh numbers for various spectral methods (SCP, SLP, CC) as well as the present method and two sets of L and a are exhibited in Table 2 . Note that SCP, SLP, and CC imply respectively the shifted Chebyshev polynomial method, the shifted Legendre polynomial method and Chebyshev collocation method [3] . Relative errors have been shown and compared with various numerical approaches in Table 2 . We noticed that compared to the Chebyshev and Legendre Spectral (SCP and SLP) approaches, the relative errors computed by our current method are significantly smaller in magnitude. All the nine critical Rayleigh numbers computed by our proposed scheme converge on the analytical ones with fewer degrees of polynomials ( n = 10). However, the parameter R achieved by the collocation method is lowered severely in the CC method [3] . From the comparison, we conclude that our present scheme is much more efficient and competes well with the other available techniques. Graphical representations of relative errors are shown in Fig. 1 . Moreover, it is eminent that non-normality is accountable for a high spectral sensitivity. The measures of a non-normality of a complex square matrix signify as non-normality ratio ( ) and can be measured by the formula given as [3] :  In this example, the non-normality ratio for the eighth order derivative matrix attained by our proposed Galerkin MWR executing the formula given in Eq. (27) is around 0.0009. Moreover, this ratio is around 0.9 with respect to the eigenvalue and does not depend upon the degree of polynomials in the range 10 < < 40 . Fig. 2 displays the curve log (Cond( A )) versus 2 < < 9 for the degree of polynomials n = 20 and n = 50 of Example 2 . It has been observed that for n = 50, the condition numbers of the eigenvalues for different values of wave numbers using twenty polynomials vary between ( 10 5 ) to ( 10 8 ) and using fifty polynomials vary between ( 10 35 ) to ( 10 37 ) and are growing for the increasing . For n = 20, the first curve moderately decreases as the wave numbers increase. It is noticed that for = 50 , ℎ second curve is almost horizontal for 6 ≤ ≤ 9 . This shows that for larger n and for various distinct values of a , the condition number ( | . . |) of the Rayleigh numbers does not worsen, indicating that the existing technique is stable. Fig. 1 demonstrates that the suggested approach's errors for the nine critical Rayleigh numbers are of order 10 −8 and, ( = 10 ) , whereas the smallest errors reported by [3] by exploiting those of Chebyshev and Legendre Spectral (SCP and SLP) methods is of order 10 −6 . As a result, the current technique achieves precisely accurate lower Rayleigh numbers.

Example 3.
We consider the eighth order ordinary differential equation worked out in [1] (28a) Fig. 2. The conditioning of matrix A with respect to the wave number. It is necessary to enforce the corresponding free-free boundary conditions given by ( 2 ) where 1 = 4 2 , 2 = 6 4 − 2 ( 2 1 + 1 ) + , 3 = 4 6 − 2 2 2 ( 2 1 + 1 ) We consider the eighth-order differential eigenvalue problem (1) in which = is purely imaginary and takes the form as given in ( 28a ). Numerical results are obtained using n = 10, 1 = 0 . 025 for the first six critical values of problem 3 are displayed in Table 3 .
Following [ 1 , 2 ], the minimum values of and the associated values of and ( = ) are computed for five distinct values of the Taylor number . For given values of , and , each eigenvalue obtained by our present method is real whereas the same are complex as demonstrated in [ 1 , 2 ]. The results are depicted in Table 3 in comparison with other numerical techniques show that the present method is compatible and much more efficient. The minimum eigenvalue is produced corresponding to the wave number which are demonstrated in [2] and [1] , respectively. It is observed that in accordance with the result given in [2] along the various over stable solutions → 1 2 as → ∞. The explicit form of the asymptotic behaviors of the Rayleigh number and the wave number of the associated disturbance for the onset of over stability may be noted; they are illustrated in detail in [2] . Note that asymptotic relations are valid only for 2 → ∞, not for → ∞, a distinction which is important when → 0 . In Fig. 3 , we illustrate the ( , )− relations for a spinning bottom layer of fluid heated below. It is noticed that for each value of < 0 . 667 , the instability sets in as ordinary cellular convection for less than a certain ( ) while it sets in as over stability for > .The values of ( ) can be determined using the conditions given for critical and minimum Rayleigh numbers discussed in [2] .

Example 4.
We examine the Sturm-Liouville problem of tenth order [14] . To which the various curves are shown   Fig. 3. ( , ) -relation for different Prandtl number for convection and over stability prediction.

Table 4
The absolute errors for first six eigenvalues of Example 4 . The first six eigenvalues utilizing Galerkin MWR exploiting for Eqs. (29a) and ( 29b ) using Bernstein polynomials along with their absolute error with ADM [14] are displayed in Table 4 . Table 4 reveals that the absolute errors | − | for all six eigenvalues are exceedingly small i.e., 1.0e + 00. This shows that our proposed method agrees well with ADM [14] .
The corresponding free-free boundary conditions are given as follows: According to [2] , here we note that the non-dimensional number plays for problem involving magnetic fields the same role which the Taylor number plays for problems involving rotation. From Eq. (2) it follows that, instability appears first for the lowest mode = 1 and the equivalent formula for is [2] where the relationship between , 1 and , 1 are can be evaluated directly as a function of (in accordance with Eq. (31) and locate the minimum numerically [2] . The eigenvalue problem in Eqs. (30a) and ( 30b ) is solved using the formulation illustrated in this section. To compare the computed results  We have noticed from the results in Table 5 The dependence of the Rayleigh number for the onset of instability on the wave number when 1 = 10 5 and 1 = 40 , 80 , 100 and 200 . The abscissa is related to the wave number as given in ( 31 a). We observe the occurrence of two minima and their relative disposition displayed in Fig. 5 . In other words, if we start with an initial situation in which 1 = 10 5 and no magnetic field exists; gradually strengthen the magnetic field. Example 6. We studied the following twelfth order eigenvalue problem [1] ( 12 − 1 10 + 2 8 − 3 6 + 4 4 − 5 2 + 6 ) ( ) + ( − 7 10 + 8 8 − 9 6 + 10 4 − 11 2 + 12 ) ( ) The corresponding free-free boundary conditions illustrated in [1] : The coefficients ( = 1 , 2 , 3   Here 1 = 2 , 1 = 4 and = 1 2 . The numerical results for using Eqs. (32a) and ( 32b ) are computed using Galerkin MWR for increasing Taylor numbers and the same magnetic field strength are displayed in Tables 6(a) , 6(b) , and 6(c) . Similar outcomes are also described in [1] . Tables 6(a) , 6(b) , and 6(c) show that our current method produces extremely similar results, except for Table 6(c) , where the results are slightly deviated from the reported results. We conclude that for higher Taylor numbers, the computed results become slightly less accurate than those of [1] .

Conclusions
In this study, we attempt a novel numerical method for computing approximate eigenvalues/critical numbers of regular eighth, tenth and twelfth order eigenvalue problems. All the derivative boundary conditions have been incorporated directly in the residual equation without reducing the order of the derivatives, with the help of orthogonal Bernstein polynomials as bases.
Our results by Galerkin MWR for 8th order problems are demonstrated in Tables 1 and 2 . In terms of computational efficiency, the numerical results reveal that the presented approach performs well with the Magnus expansion method. The comparability of the calculated Rayleigh numbers with the analytical solutions reveals a good estimation as it guarantees compatibility with the other existing methods. From Tables 3 , 4 , Table 5 and Table 6, it has been noticed that, estimated critical Rayleigh numbers corresponding to the fixed wave numbers are much closer to [ 1 , 2 ]. Thus, the existing method shows a good performance and achieves high precision when applied to numerous problems presented in this study. Although the estimated critical values of A differ slightly from those in Ref. [1] , Tables for twelfth order eigenvalue problems reveal that results produced are lower, implying that over stability is the cause of instability. This method produces a sparse co-efficient matrix with a symmetric banded matrix, which reduces computation cost and results in a well-conditioned matrix. Hence implementing the boundary conditions is much simpler and easier in these problems. The linear stability of the electro-hydro dynamic (EHD) problem was efficiently solved by exploiting Galerkin MWR. Our numerical calculations confirms that Galerkin matrix behaves better than both the Chebyshev Tau and collocation matrices. The current scheme produces a more normal discretization matrix as it produces reduced non-normality ratio. Also the resulting matrix in the latter case is also symmetric and agrees well with the other published results. The stability or instability of hydrodynamic and hydro-magnetic systems can be comprehended by a set of non-dimensional parameters. A comparisons of critical Rayleigh numbers calculated by the proposed method and those of other existing techniques are made. The obtained results confirm that the proposed method possesses much more accurate results, is more expedient and competent for higher order eigenvalue problems. Despite being more precise and effective with constant coefficients than those of the problems with variable coefficients regardless of order, the current method has a drawback in that it is less accurate and efficient with variable coefficients.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article.

Acknowledgement
The second author was partially supported by a University Grants Commission, University of Dhaka funding for the academic year 2021-2022. dimensionless vertical co-ordinate, = ( ) ∶ vertical co-ordinate.